##############################################################
######### DIFF-IN-DIFFS ANALYSIS WITH TERM LIMITS ############
######## STRICKLAND AND CROSSON, K STREET ON MAIN ############
##############################################################


setwd("C:/Users/jcrosson/Dropbox/Candidacy/K Street on Main")

require(tidyverse)
require(foreign)



### Read in Data

data <- read.csv("K_Street_on_Main_dataset_50_states.csv")
data$density00 <- data$density*100
data$biennium <- ifelse(data$year %in% c(1986), 1,
                        ifelse(data$year %in% c(1987, 1988), 2,
                               ifelse(data$year %in% c(1989, 1990), 3,
                                      ifelse(data$year %in% c(1991, 1992), 4,
                                             ifelse(data$year %in% c(1993, 1994), 5,
                                                    ifelse(data$year %in% c(1995, 1996), 6,
                                                           ifelse(data$year %in% c(1997, 1998), 7,
                                                                  ifelse(data$year %in% c(1999, 2000), 8,
                                                                         ifelse(data$year %in% c(2001, 2002), 9,
                                                                                ifelse(data$year %in% c(2003, 2004), 10,
                                                                                       ifelse(data$year %in% c(2005, 2006), 11,
                                                                                              ifelse(data$year %in% c(2007, 2008), 12,
                                                                                                     ifelse(data$year %in% c(2009, 2010), 13,
                                                                                                            ifelse(data$year %in% c(2011, 2012), 14,
                                                                                                                   ifelse(data$year %in% c(2013, 2014), 15,
                                                                                                                          16)))))))))))))))

### Table 2 Regressions ###

## Models 

require("multiwayvcov")
data$bowen000 <- data$bowen
data$genexp000000 <- data$genexp / 1000000

# Robust SE Function from James Martherus

robustse.f <- function(model, cluster, df_correction) {
  require(sandwich)
  require(lmtest)
  require(multiwayvcov)
  if(missing(cluster)) {
    name <- deparse(substitute(model))
    modelname <- paste(name,"rob",sep=".")
    model$se <- coeftest(model, vcov=vcovHC(model,"HC1"))[,2]
    model$vcovHC <- vcovHC(model,"HC1")
    assign(modelname,model,envir = .GlobalEnv)
    coeftest(model, vcov=vcovHC(model,"HC1"))
  } else {
    name <- deparse(substitute(model))
    modelname <- paste(name,"rob",sep=".")
    vcovCL <- cluster.vcov(model, cluster, df_correction = df_correction)
    model$vcovCL <- vcovCL
    modelname <- paste(name,"clustrob",sep=".")
    model$se <- coeftest(model, vcovCL)[,2]
    assign(modelname,model,envir = .GlobalEnv)
    coeftest(model, vcovCL)
  }
}

require("multiwayvcov")

model_lm_twfe <- lm(density00 ~ termlimit+ nosession + firmreport + employees + noexpire  + factor(fips) + factor(year), data = data)

model_lm_twfe_controls <- lm(density00 ~ termlimit + definitions  + prohibitions + reports + defpro + defrep  + legislators+ nosession + firmreport + employees + noexpire +foldedranney6yr+ initiativestate + bowen000 + factor(fips) + factor(year), data = data)

lt_nested <- lm(density00 ~ termlimit + nosession + firmreport + employees + noexpire + factor(fips)*year, data = data)
lt_nested.rob <- robustse.f(lt_nested)

lt_control <- lm(density00 ~ termlimit + legislators+ firmreport +foldedranney6yr + initiativestate + bowen000 +
                   nosession  + definitions + prohibitions + reports + defpro + defrep + employees + noexpire + factor(fips)*year, data = data)
lt_control.rob <-  robustse.f(lt_control)

## Residuals model (Model 9)

data <- subset(data, !(year > 2002 & abbr == "OR") ) # Remove states that switched back and fourth

base <- lm(density00 ~ legislators+ firmreport +foldedranney6yr + initiativestate + bowen000 +definitions + prohibitions + reports + defpro + defrep + employees + noexpire +
             nosession  +factor(fips) + factor(year), data = data)
predictions <- cbind(predict(base), base$model)
names(predictions)[1] <- "predicted_density"

predictions <- cbind(predictions,data.frame(residuals(base)))

# Add "ever term limit" variable

out <- NA
for(i in unique(data$abbr)){
  state <- i
  termlimit <- subset(data, abbr == i)$termlimit
  term_yes <- ifelse(sum(termlimit) > 0,1,0)
  out <- rbind(out, c(state,term_yes))
}
out <- data.frame(out)
names(out) <- c("abbr", "ever_term_limits")

data <- left_join(data, out)
data$ever_term_limits <- as.numeric(as.character(data$ever_term_limits))

cw <- data %>% dplyr::select(fips, year, ever_term_limits, termlimit)
names(predictions)[16:17] <- c("fips", "year")
predictions$fips <- as.numeric(as.character(predictions$fips))
predictions$year <- as.numeric(as.character(predictions$year))
predictions <- left_join(predictions, cw, by = c("fips", "year"))

second_stage <- lm(residuals.base. ~ factor(termlimit), data = subset(predictions, ever_term_limits == 1))
summary(second_stage)

## Write out table

require("stargazer")
require("margins")
summary(lt_nested.rob)
writeLines(capture.output(stargazer(model_lm_twfe, model_lm_twfe_controls,lt_nested.rob, lt_control.rob, second_stage,
          star.cutoffs = c(0.05, 0.01, 0.001))), "Table2.tex")


### Figure 4 ###

# Note that this code changes to a custom font for WINDOWS ONLY
# If you are not operating on a Windows machine, simply skip these formatting steps

require("ggplot2")
require("ggthemes")
require("extrafont")
loadfonts(device = "win")
require("ggrepel")
require("sysfonts")
require("showtext")
font_add_google(name = "Barlow", family = "barlow")
showtext_auto()
require(tidyverse)
require("gtools")
require("RColorBrewer")

theme_safeskies <- theme_bw()+
  theme(
    plot.margin = unit(c(.5, .5, .5, .5), units="line"),
    axis.text.x = element_text(size = 18, color="#934f3f", family = "Gill Sans MT"),
    axis.text.y = element_text(size = 15, color="#934f3f", family = "Gill Sans MT"),
    axis.title.y = element_text(size = 18, color="#934f3f", family = "Gill Sans MT"),
    plot.title = element_text(size=20, hjust=-0.03, vjust=5, margin=margin(b=15, t=9), family="Corbel", face = "bold"),
    plot.subtitle = element_text(size=9, hjust=-0.05, vjust=6, margin=margin(b=15, t=-4), family="Corbel"),
    panel.border = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.grid.major.y = element_line(size=.1, color="grey66"),
    panel.grid.minor.y = element_blank()
  )

## Data prep for Figure 4

data$density <- log(data$density00)# Logged for visualization

model_lm_twfe_controls_clients <- lm(density ~ termlimit  + legislators+ firmreport + nosession +foldedranney6yr + initiativestate + bowen000 +
                                       nosession + clients + factor(fips) + factor(year), data = data)
newdata <- model_lm_twfe_controls_clients$model
newdata$termlimit <- ifelse(newdata$termlimit == 0, 1, 0)
names(newdata)[10:11] <- c("fips", "year")
predictions <- data.frame(predict(model_lm_twfe_controls_clients, newdata = newdata))
newdata <- cbind(newdata, predictions)
names(newdata)[12] <- "alt_prediction"

out <- NA
for(i in unique(data$fips)){
  state <- i
  abbr <- as.character(subset(data, fips == i)$abbr[1])
  termlimit <- subset(data, fips == i)$termlimit
  term_yes <- ifelse(sum(termlimit) > 0,1,0)
  out <- rbind(out, c(state,abbr,term_yes))
}
out <- data.frame(out)
names(out) <- c("fips", "abbr","ever_term_limits")

out$fips <- as.numeric(as.character(out$fips))
newdata$fips <- as.numeric(as.character(newdata$fips))
newdata <- left_join(newdata, out)
parse_states <- subset(newdata, year == 2011 & is.na(alt_prediction) == F)$fips


## Code for Figure 4:

p <- ggplot(subset(data, year == 1989), aes(factor(year), exp(density)))
pdf("Figure4.pdf", 10, 8)
p +
  geom_boxplot(color = "gray", size = 1, outlier.shape = NA, width = .15, position = position_nudge(x = -.2)) + 
  stat_boxplot(geom = 'errorbar', size = 1, color = "grey", width = .15, position = position_nudge(x = -.2)) +
  geom_boxplot(data = subset(data, year == 2011), aes(factor(year), exp(density)), color = "gray", size = 1, outlier.shape = NA, width = .15, position = position_nudge(x = .2)) + 
  stat_boxplot(data = subset(data, year == 2011), aes(factor(year), exp(density)), geom = 'errorbar', size = 1, color = "grey", width = .15, position = position_nudge(x = .2)) +
  geom_point(data = subset(newdata, ever_term_limits == 0& (year == 1989 | year == 2011)), aes(year, exp(density)),fill = "darkgray", colour= "darkgray", alpha = .4, position = position_jitter(width = .08), size = 3, shape = 16) + 
  geom_point(data = subset(newdata, ever_term_limits == 0 & year == 2011),aes(year, exp(alt_prediction)), alpha = .8, fill = "black", alpha = .4, position = position_jitter(width = .08), size = 3, shape = 16) +
  geom_segment(data = subset(newdata, ever_term_limits == 0& year == 1989 & fips %in% parse_states), aes(x = factor(1989), y = exp(subset(newdata, ever_term_limits == 0& year == 1989 & fips %in% parse_states)$density),xend = factor(2011),yend = exp(subset(newdata, ever_term_limits == 0& year == 2011)$density)), colour = "darkgray", alpha = .4) +
  geom_segment(data = subset(newdata, ever_term_limits == 0& year == 1989 & fips %in% parse_states), aes(x = factor(1989), y = exp(subset(newdata, ever_term_limits == 0& year == 1989 & fips %in% parse_states)$density),xend = factor(2011),yend = exp(subset(newdata, ever_term_limits == 0& year == 2011)$alt_prediction)), colour = "black", alpha = .5) +
  theme_safeskies +
  xlab("") + ylab("Lobby Network density") 
dev.off()


################################################
##### APPENDIX MATERIAL ########################
################################################

### Figure A5 ###

require("ggjoy")

p <- ggplot(data, aes(density00, y = factor(year)))
pdf("density_year_joyplot.pdf", 8, 15)
p + geom_density_ridges2() + xlim(c(0,6)) +
  theme_safeskies +
  ggtitle("") + xlab("Network Density (x100)") + ylab("Year") +
  theme(plot.title = element_text(face="bold"), axis.text=element_text(size=11), axis.title = element_text(size = 11),
        panel.background = element_rect(fill = "transparent", colour = NA), plot.background = element_rect(fill = "transparent", colour = NA),
        legend.title = element_blank()) 
dev.off()


### First Parallel Trends Graph ###

## Figure out when each state moved in/out of treatment

terms <- subset(data, ever_term_limits == 1)

out2 <- NA
for(i in unique(terms$abbr)){
  last_year <- subset(terms, abbr == i & termlimit == 0)$year[which.max(subset(terms, abbr == i & termlimit == 0)$year)]
  out2 <- rbind(out2, c(i, last_year))
}
out2 <- data.frame(out2)
names(out2) <- c("abbr", "year_treated")
out2$year_treated <- as.numeric(as.character(out2$year_treated)) + 1

data <- left_join(data, out2)
data$year_treated <- as.numeric(as.character(data$year_treated))
data$year_treated[is.na(data$year_treated) == T] <- 0
data$biennium_treated <- ifelse(data$year_treated %in% c(1986), 1,
                                ifelse(data$year_treated %in% c(1987, 1988), 2,
                                       ifelse(data$year_treated %in% c(1989, 1990), 3,
                                              ifelse(data$year_treated %in% c(1991, 1992), 4,
                                                     ifelse(data$year_treated %in% c(1993, 1994), 5,
                                                            ifelse(data$year_treated %in% c(1995, 1996), 6,
                                                                   ifelse(data$year_treated %in% c(1997, 1998), 7,
                                                                          ifelse(data$year_treated %in% c(1999, 2000), 8,
                                                                                 ifelse(data$year_treated %in% c(2001, 2002), 9,
                                                                                        ifelse(data$year_treated %in% c(2003, 2004), 10,
                                                                                               ifelse(data$year_treated %in% c(2005, 2006), 11,
                                                                                                      ifelse(data$year_treated %in% c(2007, 2008),12, 0))))))))))))


## first treatment

data$treat_period <- ifelse(data$year < 1996, 0, 1)

data$treat_period2 <- ifelse(data$year < 1998, 0, 1)

data$treat_period3 <- ifelse(data$year < 2000, 0, 1)

data$treat_period3.5 <- ifelse(data$year < 2001, 0, 1)

data$treat_period4 <- ifelse(data$year < 2004, 0, 1)

dat <- data %>% group_by(year_treated) %>% summarise(numstates = n_distinct(fips))

data <- data %>% dplyr::group_by(fips) %>% 
  mutate(
    density_scaled = scale(density00)
  )

data$ever_tl_text <- ifelse(data$ever_term_limits == 1, "Treated", "Not Treated")

## Graph:

p <- ggplot(data = data, aes(year, density_scaled))
pdf("parallel_trends_graph1.pdf", 10, 8)
p + geom_point(aes(color = factor(ever_tl_text)), alpha = .7, shape = 16) +
  scale_color_manual(values = c("black", "gray44")) +
  stat_smooth(data = subset(data, treat_period3 == 0), aes(color = factor(ever_tl_text), linetype = factor(ever_tl_text)), method = "lm", fill = "lightgray") +
  stat_smooth(data = subset(data, treat_period3 == 1), aes(color = factor(ever_tl_text), linetype = factor(ever_tl_text)), method = "lm", fill = "lightgray") +
  geom_vline(xintercept = 2000) +
  theme_safeskies +
  ggtitle("Examining Parallel Trends Assumption, Treatment in 2000") + xlab("") + ylab("Lobby Network Density") +
  theme(plot.title = element_text(face="bold"), axis.text=element_text(size=11), axis.title = element_text(size = 11),
        panel.background = element_rect(fill = "transparent", colour = NA), plot.background = element_rect(fill = "transparent", colour = NA),
        legend.title = element_blank())
dev.off()

### Second Parallel Trends Graph (1998)

p <- ggplot(data = data, aes(year, density_scaled))
pdf("alternative_parallel_trends.pdf", 14,10)
p + geom_point(aes(color = factor(ever_tl_text)), alpha = .7, shape = 16) +
  scale_color_manual(values = c("black", "gray44")) +
  stat_smooth(data = subset(data, treat_period2 == 0), aes(color = factor(ever_tl_text), linetype = factor(ever_tl_text)), method = "lm", fill = "lightgray") +
  stat_smooth(data = subset(data, treat_period2 == 1), aes(color = factor(ever_tl_text), linetype = factor(ever_tl_text)), method = "lm", fill = "lightgray") +
  geom_vline(xintercept = 1998) +
  theme_safeskies +
  ggtitle("Examining Parallel Trends Assumption, Treatment in 1998") + xlab("") + ylab("Lobby Network Density") +
  theme(plot.title = element_text(face="bold"), axis.text=element_text(size=11), axis.title = element_text(size = 11),
        panel.background = element_rect(fill = "transparent", colour = NA), plot.background = element_rect(fill = "transparent", colour = NA),
        legend.title = element_blank()) 
dev.off()


### Table A5 ###

pred <- lm(density00 ~ nosession + firmreport + noexpire + employees, data = data)
predictions <- predict(pred)
merge_pred <- c(names(predictions))
predictions <- data.frame(cbind(merge_pred, data.frame(predictions)))

data$merge_pred <- as.character(1:nrow(data))
data <- left_join(data, predictions)
data$density00_adjusted <- data$density00 - data$predictions

require("did")

data_b <- data %>% group_by(fips, biennium) %>% summarise(biennium_treated = mean(biennium_treated), legislators = mean(legislators),density00 = mean(density00),
                                                          firmreport = mean(firmreport),density00_adjusted = mean(density00_adjusted), foldedranney6yr = mean(foldedranney6yr), 
                                                          initiativestate = mean(initiativestate),bowen000 = mean(bowen000), definitions = mean(definitions), prohibitions = mean(prohibitions),
                                                          reports = mean(reports), employees = mean(employees), noexpire = mean(noexpire))

set.seed(2452)
attgt_b <- att_gt(yname = "density00_adjusted", tname = "biennium", idname = "fips", gname = "biennium_treated",
                  data = data_b, allow_unbalanced_panel = T)
effects_b <- aggte(attgt_b, type = "dynamic", na.rm = T)
tablea5 <- summary(effects_b)


### Table A6 ###

data <- read.csv("K_Street_on_Main_dataset_50_states.csv")
data$density00 <- data$density*100

out <- NA
for(i in unique(data$abbr)){
  state <- i
  termlimit <- subset(data, abbr == i)$termlimit
  term_yes <- ifelse(sum(termlimit) > 0,1,0)
  out <- rbind(out, c(state,term_yes))
}
out <- data.frame(out)
names(out) <- c("abbr", "ever_term_limits")

data <- left_join(data, out)
data$ever_term_limits <- as.numeric(as.character(data$ever_term_limits))


pts <- lm(density00 ~ ever_term_limits*factor(year) + factor(fips), data = data)

writeLines(capture.output(stargazer(pts)), "TableA6.tex")








